! VSemiletov, 08.2011
subroutine OutputTransducerDataInit()
use Data_cfoam, only: processor, ie, nump,rank, ierr, &
                      U,p,T
use DataTime, only: ttime
use Data_Mesh, only: nCells, nPoints, meshCellCenters, meshPoints, meshCellPoints, meshCells,meshFaces
use Data_Transducer
implicit none
include 'mpif.h'
    
    integer :: iProc, iTrans, fout
    real(8),allocatable :: Send(:)
    integer :: i,j, nFL
    character(280) :: fileName

    integer :: OpenFile, ionPoints(1:nump),ionCells(1:nump)
    
    integer,allocatable :: meshCells2Face(:,:)
    integer :: gfpos(1:nTransducers*nump),fpos(1:nTransducers), &
               dpextent, iextent, secdata, secdata2, offr2(0:4), &
               secdata3(4,Transducer(1)%nFacesOut)
    integer(kind=MPI_OFFSET_KIND) :: fdispl
!    double precision :: rsecdata(1:3,1:nPoints)
    
    ierr=0

!    call MPI_GET_ADDRESS(firstdata(1)%item1, offr2(0), ierr) 
!    write(*,*) "rank 2", offr2(0), offr2(1),offr2(2),offr2(3),offr2(4)  

!   addresses of individual components of the datatype are not contigous in
!   memory due to allocatable memory, therefore not possible to use
!   MPI_TYPE_STRUCT with a single IO operation - need individual ones to 
!   different memory locations

! calculate and  setup offsets
!    call MPI_TYPE_EXTENT deprecated 
    call MPI_TYPE_SIZE(MPI_INTEGER, iextent, ierr)
    call MPI_TYPE_SIZE(MPI_DOUBLE_PRECISION, dpextent, ierr)

!   sending arrays definition
!   calculate offsets for MPI-IO
    fpos=0
    do i = 1,nTransducers
      fpos(i)=fpos(i)+ (  (3*nPoints*dpextent) + (8*Transducer(i)%nCellsOut*iextent) & 
                         +(4*Transducer(i)%nFacesOut*iextent) + 3*iextent )
    end do

    call MPI_Allgather( fpos, nTransducers, MPI_INT, gfpos, &
                        nTransducers, MPI_INT, MPI_COMM_WORLD,ierr)
    
    call MPI_Barrier(MPI_COMM_WORLD,ierr)

    fdispl=0

    if (rank>0) then
    do j = 0,rank
      do i = 1, nTransducers
        fdispl = fdispl + gfpos(nTransducers*(j-1)+i)
      end do
    end do
    end if

!   end of calculation for the offsets

    do iTrans = 1,nTransducers
        if(Transducer(iTrans)%AllProcsWriteOwnData)then
            ierr=0
            call DefineNextFileName(Transducer(iTrans)%fileOutput,Transducer(iTrans)%nameLength,.false.)
            call MPI_FILE_OPEN(MPI_COMM_WORLD,Transducer(iTrans)%fileOutput(1:Transducer(iTrans)%nameLength), &
                               MPI_MODE_WRONLY + MPI_MODE_CREATE, MPI_INFO_NULL, fout, ierr) 

!           the following lines open up a new file for each process
            Transducer(iTrans)%nameLength=Transducer(iTrans)%nameLength - 7
            call DefineNextFileName(Transducer(iTrans)%fileOutput,Transducer(iTrans)%nameLength,.true.)
            open(1,file=Transducer(iTrans)%fileOutput(1:Transducer(iTrans)%nameLength),form='unformatted',IOSTAT=OpenFile)
            if (OpenFile.ne.0) then
                write(*,*) '@@ Fatal Error: File open error ', Transducer(iTrans)%fileOutput(1:Transducer(iTrans)%nameLength)
                STOP '@@ Fatal Error!'
            end if
            close(1)
!           end of lines for opening up a new file for each process

            call MPI_FILE_CLOSE(fout, ierr) 
            Transducer(iTrans)%nCellsOutGlobal = 0
            call MPI_ALLREDUCE(Transducer(iTrans)%nCellsOut,Transducer(iTrans)%nCellsOutGlobal,1,MPI_INTEGER, &
                               MPI_SUM,MPI_COMM_WORLD,ierr)
        
            if(rank.eq.0)then
                write(*,*) 'Trans ',iTrans,' nCellsGlobal=', Transducer(iTrans)%nCellsOutGlobal
            end if
        
            cycle
        end if

!        write(*,*) '@@@ nTransducers=',iTrans, rank
        if(rank.ne.0)then
            allocate(Transducer(iTrans)%recvCounts(1),Transducer(iTrans)%rDispls(1))
            allocate(Transducer(iTrans)%recvCountsFace(1),Transducer(iTrans)%rDisplsFace(1))
        else
            allocate(Transducer(iTrans)%recvCounts(nump),Transducer(iTrans)%rDispls(nump+1))
            allocate(Transducer(iTrans)%recvCountsFace(nump),Transducer(iTrans)%rDisplsFace(nump+1))
        end if
    
        call DefSendingArrays(Transducer(iTrans)%nCellsOutGlobal, &
                              Transducer(iTrans)%nCellsOut, &
                              Transducer(iTrans)%recvCounts, &
                              Transducer(iTrans)%rDispls) !, &
                              !Transducer(iTrans)%recv)

        allocate(Transducer(iTrans)%recv(Transducer(iTrans)%nCellsOutGlobal))

        call DefSendingArrays(Transducer(iTrans)%nFacesOutGlobal, &
                              Transducer(iTrans)%nFacesOut, &
                              Transducer(iTrans)%recvCountsFace, &
                              Transducer(iTrans)%rDisplsFace) !, &
                              !Transducer(iTrans)%recvFace)

        allocate(Transducer(iTrans)%recvFace(Transducer(iTrans)%nFacesOutGlobal))

        if(rank.eq.0)then
            write(*,*) 'Trans nCellsGlobal=', Transducer(iTrans)%nCellsOutGlobal
            write(*,*) 'Trans nFacesGlobal=', Transducer(iTrans)%nFacesOutGlobal
        end if
        
        call DefineNextFileName(Transducer(iTrans)%fileOutput,Transducer(iTrans)%nameLength,.false.)
        
        open(1,file=Transducer(iTrans)%fileOutput(1:Transducer(iTrans)%nameLength),form='unformatted',IOSTAT=OpenFile)
        if (OpenFile.ne.0) then
            write(*,*) '@@ Fatal Error: File open error ', Transducer(iTrans)%fileOutput(1:Transducer(iTrans)%nameLength)
            STOP '@@ Fatal Error!'
        end if
        close(1)
    end do

    do iTrans = 1,nTransducers

        ! Mesh output for Slice:    
        if(Transducer(iTrans)%SliceOut)then
            
            if(Transducer(iTrans)%AllProcsWriteOwnData)then
                nFL = Transducer(iTrans)%nameLength - 7
            else
                nFL = Transducer(iTrans)%nameLength
            end if
            fileName(1:nFL) = Transducer(iTrans)%fileOutput(1:nFL)
            fileName(nFL+1:nFL+4) = 'mesh'
            nFL = nFL+4

            if(Transducer(iTrans)%AllProcsWriteOwnData)then
                open(1,file=fileName(1:nFL),form='unformatted')
                ierr=0
                nFL = Transducer(iTrans)%nameLength - 7
                fileName(1:nFL) = Transducer(iTrans)%fileOutput(1:nFL)
                fileName(nFL+1:nFL+4) = 'mesh'
                nFL = nFL+4
                i=LEN(processor(1:ie))
!               the filename will need to be changed if not writing to
!               individual files as well
                call MPI_FILE_OPEN(MPI_COMM_WORLD,fileName(1+i:nFL), &
                                   MPI_MODE_WRONLY + MPI_MODE_CREATE, MPI_INFO_NULL, fout, ierr)
                call MPI_FILE_SEEK(fout,fdispl,MPI_SEEK_SET,ierr)

            else
    	        open(1,file=processor(1:ie)//fileName(1:nFL),form='unformatted')
    	    end if
	        write(1) nPoints,Transducer(iTrans)%nCellsOut, &
	                        ((meshPoints(j,i),j=1,3),i=1,nPoints), &
                            ((meshCellPoints(j,Transducer(iTrans)%CellsOut(i)),j=1,8), &
                            i=1,Transducer(iTrans)%nCellsOut)
            write(1) Transducer(iTrans)%nFacesOut, &
!	                        ((meshCells(j,Transducer(iTrans)%CellsOut(i)),j=1,6),i=1,Transducer(iTrans)%nCellsOut), &
	                        ((meshFaces(j,Transducer(iTrans)%FacesOut(i)),j=1,4),i=1,Transducer(iTrans)%nFacesOut)
            close(1)

            call MPI_FILE_WRITE_ALL(fout, nPoints, 1, MPI_INTEGER, &
                                MPI_STATUS_IGNORE, ierr)
            call MPI_FILE_WRITE_ALL(fout, Transducer(iTrans)%nCellsOut, 1, MPI_INTEGER, &
                                MPI_STATUS_IGNORE, ierr)
            call MPI_FILE_WRITE_ALL(fout, meshPoints(1,1),3*nPoints,MPI_DOUBLE_PRECISION, &
                                MPI_STATUS_IGNORE, ierr) 
            call MPI_FILE_WRITE_ALL(fout, meshCellPoints(1,1),8*Transducer(iTrans)%nCellsOut, &
                                    MPI_INTEGER, MPI_STATUS_IGNORE, ierr)
            call MPI_FILE_WRITE_ALL(fout,Transducer(iTrans)%nFacesOut, 1, &
                                    MPI_INTEGER, MPI_STATUS_IGNORE, ierr)
            call MPI_FILE_WRITE_ALL(fout, meshFaces(1,1), 4*Transducer(iTrans)%nFacesOut, &
                                    MPI_INTEGER, MPI_STATUS_IGNORE, ierr)

            call MPI_FILE_CLOSE(fout, ierr)
 
!    i=LEN(processor(1:ie))

!    call MPI_FILE_OPEN(MPI_COMM_WORLD,fileName(1+i:nFL), &
!                       MPI_MODE_RDONLY, MPI_INFO_NULL, fout, ierr)
!    call MPI_FILE_SEEK(fout,fdispl,MPI_SEEK_SET,ierr)

!    call MPI_FILE_READ(fout, secdata, 1, MPI_INTEGER, &
!                                MPI_STATUS_IGNORE, ierr)

!    call MPI_FILE_READ(fout, secdata2, 1, MPI_INTEGER, &
!                                MPI_STATUS_IGNORE, ierr)
!    call MPI_FILE_READ(fout, rsecdata, 3*nPoints,MPI_DOUBLE_PRECISION, &
!                                MPI_STATUS_IGNORE, ierr)

!    call MPI_FILE_SEEK(fout,fdispl+(dpextent*3*nPoints)+(iextent*8*Transducer(iTrans)%nCellsOut) &
!                         + (iextent*3) , MPI_SEEK_SET,ierr)

!     call MPI_FILE_READ(fout, secdata3, 4*Transducer(iTrans)%nFacesOut, MPI_INTEGER, &
!                                MPI_STATUS_IGNORE, ierr)

!    call MPI_FILE_CLOSE(fout, ierr)

!            write(*,*) "Mesh files",fileName(1:nFL) 
            do j = 1,5
                do i = 1,Transducer(iTrans)%nCellsOut
                    Transducer(iTrans)%Vars(i,j) = 0.0d0
                end do
            end do
        end if

        ! Head of Table output:

        ! Format for output-file:
        ! The start of output file for Transduser
        ! x    ........................
        ! y    ........................
        ! z    ........................
        ! Time Pres Uvel Vvel Wvel Temp
        ! .... .... .... .... .... ....
        if(Transducer(iTrans)%TableOut)then

            allocate(Send(Transducer(iTrans)%nCellsOut))

            if(allocated(Transducer(iTrans)%recv).eqv..false.)then
                allocate(Transducer(iTrans)%recv(Transducer(iTrans)%nCellsOut))
            end if
            
            if(rank.eq.0)then
                open(1,file=Transducer(iTrans)%fileOutput) !,access='append')
                
                write(1,'(a12)',ADVANCE = 'NO') '@ ' ! \)',ADVANCE = 'NO') '@ '
                write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                do j = 1,Transducer(iTrans)%nVars
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(i15)',ADVANCE = 'NO') i ! \)',ADVANCE = 'NO') i
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end do
                write(1,'(a1)') ' '
!                write(1,'(a1/)') ' '
            end if
            
            ! write x
            do i = 1,Transducer(iTrans)%nCellsOut
                Send(i) = -meshCellCenters(1,Transducer(iTrans)%CellsOut(i))
                if(rank.eq.0)then
                    Transducer(iTrans)%recv(i) = -meshCellCenters(1,Transducer(iTrans)%CellsOut(i))
                end if
            end do
            
            if(nump.gt.1)then
                call MPI_Gatherv(Send,Transducer(iTrans)%nCellsOut,MPI_DOUBLE_PRECISION,Transducer(iTrans)%recv, &
                             Transducer(iTrans)%recvCounts,Transducer(iTrans)%rDispls, &
                             MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            end if
            
            if(rank.eq.0)then
                write(1,'(a12)',ADVANCE = 'NO') 'x ' ! \)',ADVANCE = 'NO') 'x '
                write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                do j = 1,Transducer(iTrans)%nVars
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(f15.5)',ADVANCE = 'NO') Transducer(iTrans)%recv(i) ! \)',ADVANCE = 'NO') Transducer(iTrans)%recv(i)
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end do
                write(1,'(a1)') ' '
!                write(1,'(a1/)') ' '
            end if
            
             ! write y
            do i = 1,Transducer(iTrans)%nCellsOut
                Send(i) = -meshCellCenters(2,Transducer(iTrans)%CellsOut(i))
                if(rank.eq.0)then
                    Transducer(iTrans)%recv(i) = -meshCellCenters(2,Transducer(iTrans)%CellsOut(i))
                end if
            end do
            
            if(nump.gt.1)then
                call MPI_Gatherv(Send,Transducer(iTrans)%nCellsOut,MPI_DOUBLE_PRECISION,Transducer(iTrans)%recv, &
                             Transducer(iTrans)%recvCounts,Transducer(iTrans)%rDispls, &
                             MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            end if

            if(rank.eq.0)then
                write(1,'(a12)',ADVANCE = 'NO') 'y ' ! \)',ADVANCE = 'NO') 'y '
                write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                do j = 1,Transducer(iTrans)%nVars
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(f15.5)',ADVANCE = 'NO') Transducer(iTrans)%recv(i) ! \)',ADVANCE = 'NO') Transducer(iTrans)%recv(i)
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end do
                write(1,'(a1)') ' '
!                write(1,'(a1/)') ' '
            end if

             ! write z
            do i = 1,Transducer(iTrans)%nCellsOut
                Send(i) = -meshCellCenters(3,Transducer(iTrans)%CellsOut(i))
                if(rank.eq.0)then
                    Transducer(iTrans)%recv(i) = -meshCellCenters(3,Transducer(iTrans)%CellsOut(i))
                end if
            end do
            
            if(nump.gt.1)then
                call MPI_Gatherv(Send,Transducer(iTrans)%nCellsOut,MPI_DOUBLE_PRECISION,Transducer(iTrans)%recv, &
                             Transducer(iTrans)%recvCounts,Transducer(iTrans)%rDispls, &
                             MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
            end if

            if(rank.eq.0)then
                write(1,'(a12)',ADVANCE = 'NO') 'z ' ! \)',ADVANCE = 'NO') 'z '
                write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                do j = 1,Transducer(iTrans)%nVars
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(f15.5)',ADVANCE = 'NO') Transducer(iTrans)%recv(i) ! \)',ADVANCE = 'NO') Transducer(iTrans)%recv(i)
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end do
                write(1,'(a1)') ' '
!                write(1,'(a1/)') ' '

            ! write string with titles of vars
                write(1,'(a12)',ADVANCE = 'NO') 'Time ' ! \)',ADVANCE = 'NO') 'Time '
                write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
            
                if(Transducer(iTrans)%Pres)then
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(a15)',ADVANCE = 'NO') 'Pres' ! \)',ADVANCE = 'NO') 'Pres'
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end if
                if(Transducer(iTrans)%Uvel)then
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(a15)',ADVANCE = 'NO') 'Uvel' ! \)',ADVANCE = 'NO') 'Uvel'
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end if
                if(Transducer(iTrans)%Vvel)then
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(a15)',ADVANCE = 'NO') 'Vvel' ! \)',ADVANCE = 'NO') 'Vvel'
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end if
                if(Transducer(iTrans)%Wvel)then
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(a15)',ADVANCE = 'NO') 'Wvel' ! \)',ADVANCE = 'NO') 'Wvel'
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end if
                if(Transducer(iTrans)%Temp)then
                    do i = 1,Transducer(iTrans)%nCellsOutGlobal
                        write(1,'(a15)',ADVANCE = 'NO') 'Temp' ! \)',ADVANCE = 'NO') 'Temp'
                    end do
                    write(1,'(a3)',ADVANCE = 'NO') ' | ' ! \)',ADVANCE = 'NO') ' | '
                end if
                write(1,'(a1)') ' '
!                write(1,'(a1/)') ' '
            end if
            
            deallocate(Send)

            close(1)
        end if
    end do

!    deallocate(Send)

end subroutine OutputTransducerDataInit

subroutine DefSendingArrays(nCellsOutGlobal,nCellsOut,recvCounts,rDispls) ! ,recv)
use Data_cfoam, only: processor, ie, nump,rank, ierr
implicit none
include 'mpif.h'

    integer,intent(in) :: nCellsOut
    integer,intent(out) :: recvCounts(nump), rDispls(nump+1)
    !real(8),intent(out) :: recv(nump)
!    integer,allocatable,intent(out) :: recvCounts(:), rDispls(:)
!    real(8),allocatable,intent(out) :: recv(:)
    integer,intent(out) :: nCellsOutGlobal

!    real(8),allocatable :: recv(:)
    integer :: iProc

!    write(*,*) '@@@ DefSendingArrays', rank

    if(nump.eq.1)then
!        do iTrans = 1,nTransducers
            nCellsOutGlobal = nCellsOut
!            allocate(recvCounts(1),rDispls(1))
            recvCounts(1) = nCellsOut
            rDispls(1) = 0
!            allocate(recv(nCellsOutGlobal))
!        end do
        return
    end if

        if(rank.eq.0)then
!            write(*,*) '@@@ allocate(recvCounts(nump),rDispls(nump+1))'
!            do iTrans = 1,nTransducers
!                allocate(recvCounts(nump),rDispls(nump+1))
                recvCounts = 0
!            end do
        else
!            do iTrans = 1,nTransducers
!                allocate(recvCounts(1),rDispls(1))
                recvCounts = 0
!            end do
        end if

!        do iTrans = 1,nTransducers

!            write(*,*) 'nCellsOut =', nCellsOut, rank

            call MPI_Gather(nCellsOut,1,MPI_INTEGER,recvCounts,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

            call MPI_Barrier(MPI_COMM_WORLD,ierr)
            
            if(rank.eq.0)then
                write(*,*) 'recvCounts=',recvCounts
                rDispls(1) = 0
                do iProc = 1,nump
                    rDispls(iProc+1) = rDispls(iProc) + recvCounts(iProc)
!                    write(*,*) 'Output from first proc: recvCounts(',iProc,')=', Transducer(iTrans)%recvCounts(iProc)
!                    write(*,*) 'Output from first proc: rDispls(',iProc,')=',    Transducer(iTrans)%rDispls(iProc+1)
                end do
                
                nCellsOutGlobal = rDispls(nump+1)
                
                write(*,*) 'Trans nCellsGlobal/nFacesGlobal=', nCellsOutGlobal
                
!                allocate(recv(nCellsOutGlobal))
            end if
!        end do
!    end if

end subroutine DefSendingArrays

subroutine WriteBoundaryOut(nCellsOut,CellsOut,BoundaryOut,name, InitVal,ScaleFactor,nLayers)
use Data_cfoam, only: btype, &
                     processor, ie, ierr, Boundary2DmeshOut, nump, rank
use Data_Mesh, only: nBoundaryFields, boundaryMeshStart, boundaryMeshSize, meshFaceOwner,meshFaces, &
                     meshCellCenters, meshPoints, meshFaceCenters, meshSf
implicit none
    include 'mpif.h'

    integer,intent(in) :: nCellsOut, CellsOut(nCellsOut), BoundaryOut
    character(50),intent(in) :: name
    real(8),intent(in) :: InitVal,ScaleFactor ! for 2D mesh normal to the boundary
    integer,intent(in) :: nLayers
    
    integer :: iCell,jCell,iFace, i,j, nLength, a
    integer :: nBoundCellsOut, BoundCellsOut(nCellsOut), BoundFacesOut(4,nCellsOut)
    real(8) :: r(3)
    
    real(8),allocatable :: FaceCenters(:,:), FaceNormals(:,:), FaceCentersGlobal(:,:), FaceNormalsGlobal(:,:)
    integer :: nBoundCellsOutGlobal
    real(8) :: L,b
    
    integer :: iProc
    real(8) :: normal(3), Point(3)
    
!    integer,allocatable :: FaceNumbers(:)

    if(BoundaryOut.le.0) return
    if(BoundaryOut.gt.nBoundaryFields)then
        write(*,*) '@@ Fatal error: ', name, ' BoundaryOut.gt.nBoundaryFields: ',BoundaryOut,'>',nBoundaryFields
        stop '@@ Fatal error'
    end if

    nBoundCellsOut = 0
    BoundCellsOut(:) = 0
    
    if(Boundary2DmeshOut)then
        allocate(FaceCenters(3,nCellsOut),FaceNormals(3,nCellsOut))
!        allocate(FaceNumbers(nCells))
    end if
    
!    do j=1,nBoundaryFields
!        if (btype(j).eq.BoundaryOut) then
            do iFace=1+boundaryMeshStart(BoundaryOut),boundaryMeshStart(BoundaryOut)+boundaryMeshSize(BoundaryOut)
                jCell = meshFaceOwner(iFace)
                do iCell = 1,nCellsOut
                    if(CellsOut(iCell).eq.jCell)then
                        nBoundCellsOut = nBoundCellsOut + 1
                        BoundCellsOut(nBoundCellsOut) = jCell
                        BoundFacesOut(1:4,nBoundCellsOut) = meshFaces(1:4,iFace)

                        if(Boundary2DmeshOut)then
                            FaceCenters(1:3,nBoundCellsOut) = meshFaceCenters(1:3,iFace)
                            FaceNormals(1:3,nBoundCellsOut) = meshSf(1:3,iFace)
                        end if
                    end if
                end do
            end do
!        end if
!    end do

    nLength = 50
    do j = 1,50
        if(name(j:j).eq.' ')then
            nLength = j-1
            exit
        end if
    end do

    open(11,file=processor(1:ie)//name(1:nLength),iostat=ierr)
    if(ierr.ne.0) then
        write(*,*) '@@ Fatal error: ', processor(1:ie)//name, ' open error!'
        stop '@@ Fatal error'
    end if

    write(11,*) 'nCellsOut:'
    write(11,*) nBoundCellsOut
    write(11,*) 'CellsOut  meshCellCenters(1:3)'

    do i = 1,nBoundCellsOut
        write(11,'(i15,3f20.10)') BoundCellsOut(i), &
                                meshCellCenters(1:3,BoundCellsOut(i))
    end do

    write(11,*) 'PointsOut:'
    
    do i=1,nBoundCellsOut
        r = 0.0d0
        do a = 1,4
            r(1:3) = r(1:3) + meshPoints(1:3,BoundFacesOut(a,i))
        end do
        r(1:3) = 0.25d0 * r(1:3)
        write(11,'(4i15,3f20.10)') BoundFacesOut(1:4,i), r(1:3)
    end do

    close(11)
    
! --- start output 2D mesh normal to the boundary ------------------------------------------ !
    ! j = 1: points placed on the boundary
    ! j = 2: points placed in the distance a in normal direction to the boundary
    !     a = InitVal
    ! j = 3: points placed in the distance a+a*b in normal direction to the boundary
    !     b = ScaleFactor
    ! j = 4: a + a*b + a*b^2
    !    ....

    if(Boundary2DmeshOut)then

        if(nLayers.eq.0)then
            return
        end if

        call MPI_ALLREDUCE(nBoundCellsOut,nBoundCellsOutGlobal,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)

        do iProc = 0,nump-1

            if(rank.eq.iProc)then
                if(rank.eq.0)then
                    open(11,file=trim(name)//'2DnormalMesh',iostat=ierr)
!                    write(11,*) 'VARIABLES = "x"'
!                    write(11,*) '"y"'
!                    write(11,*) '"z"'
!                    write(11,*) '"L"'
!                    write(11,*) 'ZONE'
                    
                else
                    open(11,file=trim(name)//'2DnormalMesh',position='append',iostat=ierr)
                end if
                if(ierr.ne.0) then
                    write(*,*) '@@ Fatal error: ', trim(name)//'2DnormalMesh', ' open error!'
                    stop '@@ Fatal error'
                end if
                
                if(rank.eq.0)then
                    write(11,*) 'TITLE = "2DnormalMesh"'
                    write(11,*) 'VARIABLES = "X"'
                    write(11,*) '"Y"'
                    write(11,*) '"Z"'
                    write(11,*) '"L"'
                    write(11,*) 'ZONE'
                    write(11,*) 'STRANDID=0, SOLUTIONTIME=0'
                    write(11,*) 'I=',nLayers,' J=', nBoundCellsOutGlobal, 'K=1, ZONETYPE=Ordered'
                    write(11,*) ' DATAPACKING=POINT'
                    write(11,*) 'DT=(SINGLE SINGLE SINGLE SINGLE )'
                end if

                do i = 1,nBoundCellsOut
                    L = 0d0
                    b = ScaleFactor
                    normal(1:3) = FaceNormals(1:3,i)
                    normal(1:3) = - normal(1:3) / sqrt( sum(normal(1:3)*normal(1:3)) )
                    do j = 1,nLayers
                        Point(1:3) = FaceCenters(1:3,i) + L * normal(1:3)

                        write(11,'(4e25.12)') Point(1:3), L

                        L = L + InitVal * b
                        b = b * ScaleFactor
                    end do
                end do

                close(11)
            end if

            call MPI_Barrier(MPI_COMM_WORLD,ierr)

        end do
        
        if(rank.eq.0)then
            write(*,*) '2D normal mesh (',trim(name)//'2DnormalMesh) has been written successfully!'
        end if
        
    end if
    
! -- finish output 2D mesh normal to the boundary ------------------------------------------ !

end subroutine writeBoundaryOut
